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COMPUTATIONAL METHODS FOR AERODYNAMIC DESIGN 
USING NUMERICAL OPTIMIZATION 

by 

Martin F. Peeters 
ABSTRACT 

Five methods to increase the computational efficiency of aerodyn^Lmic 
design using numerical optimization^ by reducing the computer time required 
to perform gradient calculations^ are examined. Four of these methods have 
flaws, while one shows promise. The promising method consists of drastically 
reducing the sire of the computational domain on which aerodynaunic calcula- 
tions are made during gradient calculations. Since a gradient calculation 
requires the solution of the flow adDout an airfoil whose geometry has been 
slightly perturbed from a base airfoil, the flow about the base airfoil is 
used to determine boundary conditions on the reduced computational domain. 
This method worked well in subcritical flow, but some unresolved problems 
remain if it is used in supercritical flow. 
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LIST OF SYMBOLS 

a speed of sound (ra/sec) 

c airfoil chord (m) 

CJ drag coefficient 

Cjf lift coefficient 

Cp pressure coefficient 

(0 doublet strength 

f airfoil geometry 

G constraint function 

h shape function 

K transonic similarity parameter 

Ku curvature = j' 

M Mach number 

OBJ objective function 

r coefficient of shape function 

t airfoil thickness (ra) 

u,v perturbation velocity in x and y direction, 

respectively 

area within airfoil divided by iS 
Cartesian coordinates 
dimensional values of x,y (m) 
transonic lateral coordinate 

ratio of specific heats 
thickness/choru 

pertubation velocity potential 


Vol 
x,y 
X, y 
y 
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velocity potential 


Subscripts 


b base flow 

i,j mesh point indices in x and y directions, 

respectively 


L 

local 

condition 

Is 

lower 

surface 

us 

upper 

surface 

or> 

freestream condition 


Superscripts 

q design iteration number 

' pertubation from base flow 
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Improved methods for the design of airfoils are always a 
subject of interest in aeronautical engineering. To date the 
most successful analytical mfthods for the design of 
airfoils have relied on some form of inverse calculations. 
An inverse calculation is one in which the desired flow 
field is specified and the airfoil shape is solved for, 
which generates this flow field. Examples of the use of 
inverse methods can be found in the work of Henne (ref. 1), 
or Chin and Rizzetta (ref. 2y. 

Inverse methods can be a definite aide in the design of 
airfoils but they do have some inherent drawbacks; 1) 
Inverse methods require apriori knowledge of the desired 
pressure or velocity distibution along the airfoil, 2 ) The 
desired flow field may be impossible to realize with any 
physically realistic airfoil shape, and 3) Constraints on 
the airfoils characteristics are not easy to implement. 

An alternative approach for the design of airfoils has 
been proposed by Hicks, Murman, and Vanderpl aats (ref. 3). 
The technique is to desgn airfoils using numerical 
optimization in which an aerodynamic analysis code is 
coupled to a numerical optimization code. This method would 
allow the designer to optimize a single performance 
characteristic of the airfoil while at the same time 
constraining other performance characteristics of the 
airfoil to be within certain values prescibed by the 
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designer. 

This method for designing airfoils is very flexible. It 
gives the designer freedom to choose which performance 
characteristics to optimize and how much degradation in 
other performance characteristics is tolerable. Which 
performance characteristic is to be optimized and which 
performance characteristics are to be constraints can be 
varied;, giving the designer accurate information on design 
tradeoffs. 

Initial work with this method has shown that while it 
seems to work, it requires a considerable amount of c.p.u. 
time, limiting its usefullness. The objective of this 
present research is to explore ways in which the 
computational effciency of designing airfoils using 
numerical optimization can be increased. 

The basic concepts involved in optimization will be 
reviewed first. Thereafter modifications to the method that 
could improve the computational efficiency will oe 
discussed.. Conclusions and recomendations will then be made. 

Optimization Concepts 

Consider an airfoil in which the upper surface is defined 
by the functional relationship fi/s (x/c) and the lower 
surface is defined by the functional relations lip (x/c). 
It is desired that a certain performance characteristic of 
the airfoil is optimized. For example, assume that the drag 
at zero angle of attack is to be m’nimzed. In this case the 
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drag would be called the objective function. In doing the 
optimization, it is desired that constraints be imposed on 
other performance chr., acteristics of the airfoil. For 
example, the airfoil volume is constrained to be greater 
than a specified minimum and the lift is constrained to be 
greater than a specified minimum. 

To perform the optimization, modifications will have to 
be made to the airfoil geometry. This is accomplished by 
adding shape functions to the initial airfoil so the umer 
surface of the airfoil would be given by 

fvs ( x/c) s Ui {%/ q) * r, *h, (x/c) +...+ r^*h^(x/c) 

and similarly for the lower surface. The function h„(x/c) is 

a shape function. Figure 1 illustrates some examples of 

shape functions that have been used by Hicks and 
Vanderplaats (ref. for optimization. The r^'s determine 
the magnitude of each shape function added to the initial 
airfoil. These are the only quantities that are varied in 
the optimization process so the r„ 's are refered to as the 
design variables. 

The statement of the problem can be summarized as 

follows : 

Minimize OBJ(X) 

Subject to: G; (X) < 0 i = 1,m 
where X is a design vector 
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The following terminology is useful when discussing 

optimization. The n dimensional space spanned by the vector 
X is refered to as the design space. A constraint is said 
to be inactive if G;(X) < 0; it is said to be violated if G.- 
(X) > 0 ; it is said to be active if G; (*X) sO . Since an 
exact zero is rarely found on a computer, a more reasonable 
definition for an active constraint isjG, (X)j <£ where £ is 
a small value. This will be the definition of an active 
constraint used in this paper. A design is feasible if for 
all i Gi (X) <0. Ine minimal feasible design is said to be 
optimal . 

How the optimization procedure actually works is best 
explained by illustating a simple exampl*. Consider the 
problem : 


Minimize Cj 


Subject 

to : 



1 

o 

0 (lift constraint) 

and 



Vol^,„ 

- Vol 

< 0 (airfoil volume constaint) 


/ r, \ 


where X = 


rj and r 2 are the coefficients 

two shape functions 



The design space for this hypothetical problem is shown in 
figure 2. Contours of constant objective function ( C^ ) and 
the constraints are illustrated in the design space. Assume 
that the initial airfoil is given by 
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r, r 0 
* 0 

and that this Initial airfoil is a feasible design. 

The optimization process is an iterative procedure in 
which the following recursive relationship is used 

1 s it* ♦ X * ^ C* 1 

where q is the iteration number, vector S is the direction 
of search in the design space and ^ is a scalar defining 
the distance of trav . in the direction given by S. Each 
optimization iteration thus procedcs in two steps: First the 
vector S giving the direction of travel is found, then the 
scalar is determined. 

The procedure for determining S is somewhat different 
depending on whether any constraints are active. In the 
example of figure 2, it is assumed that no constraints are 
active initially so the determination of S proceeds as 
follows. Each design variable is separately perturbed to 
determine its effect on the objective function; thus a 
finite difference approximation to the gradient of the 
objective function is constructed as 


V08J = 

I h -K. 

> COQJ) J 


f 

1 


• 

1 — — ^ 

• 


• 

a(CSJ) 


A coaj) 


^ 1 


\ h 
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In this example there are only two design variables. Fcr 
each one an eerodynamic calculation must be performed v ith 
the design variable perturbed to determine the change in the 
objective function. 

Whith a finite difference appro::imation to the gradient 
of the objective function found, S can now be determineo. 
Different optimization schemes use different methods to 
determine S. A steepest descent method would just make S the 
negative of the gradient. Conjugate gradient methods (ref. 
5) or quasi Newton methods (ref. 6) determine S as some 
function of the gradient. In general, optimization schemes 
determine S as some function of the gradient of the 
objective function. 

With S known, Ck must now be found. is found by 
conducting a linear search in the direction of S until a 
minimum is found, or until a constraint becomes active. 
Again, different optimization codes will use different 
methods to find ci . a typical method woula be to perform 3 
evaluations of the objective function on the line defined by 
S. A quadratic fit is then made with these 3 points and the 
minimum is found. Similar ideas are used in other methods. 

Now equation 2 can be used to determine a new airfoil 
geometry. In the example given in figure 2 this would cause 
a movement in the design space from A to B. This procedure 
is repeated until either a minimum is found in the objective 
function, or a constraint becomes active. In the example of 
figure 2 the optimization procedure would move the design to 
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C where the lift constraint becomes active. 

At this point a somewhat different procedure is used to 
find the vector S. In addition to finding a finite 
difference approximation to the gradient of the objective 
function, a finite difference approximation to the gradient 
of the constraint function must be found. This is done in 
the same manner as before. Each design variable is perturbed 
separately to determine its effect and the objective and 
constraint functions. The determination of S again varies 
from program to program, but usually the optimization 
program will try to move the design closer to a minimum and 
at the same time push it away slightly from the constraint. 

Constrained and unconstrained optimization iterations are 
performed as necessary until a feasible design is found 
which minimzes the objective function (optimal design). In 
the example of figure 2, this would move the design to point 

E. There is no guarantee that this design is at an absolute 
minimum. There may be rany local minima in the design space. 
In general an improvement in the design will have been made, 
but to have more assurance in finding the absolute minimum 
the optimization procedure should be started at different 
points in the design space. In the example of figure 2 point 
E is a relative minimum. If the procedure was begun at point 

F, the absoulute minimum, point G, would have been found. 


OF POOR QuKLil'V 

Examples of Airfoil Optimization 

The concepts outlined above have been tested by numerous 
people. Numerical optimization has been used to optimize low 
speed, high lift airfoils (ref. 7) ; it has been used to 
optimize airfoils in transonic flow (ref. 8,9). The results 
obtained by Hicks, Murman ana Vanderplaats (ref. 3) in which 
airfoils in transonic flow are optimized, will be presented 
as an example demonstrating the potential that this method 
has for the design of airfoils. 

Their optimization procedure couplea together an 
aerodynamic analysis code based on the small disturbance 
transonic potential equation and CONMIN (ref. 10), a FORTRAN 
program for constrained function minimization. 

Some results from their work are presented in figure 3- 
In each case the objective was to minimze drag (the only 
drag present in this inviscid calculation is wave drag ). In 
each case the freestream Mach number was 0.8; there were 
seven design variables; the airfoils were symmetric. In the 
cases of figures 3a and 3b the only constraint was an 
airfoil volume constraint, the cases of figures 3c and 3d 
imposed a curvature constraint on the airfoil and a 
thickness/chord constraint in addition to the airfoil volume 
constraint. In every case a significant reduction in drag 
was realized. 
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Limitations of Airfoii Design Using Numerical Optimization 

The above examples indicate that improved airfoil 
designs can be realized using numerical optimization as a 
design tool. This method is, however, limited by its 
excessive appetite for c.p.u. time. The examples presented 
utilized an aerodynamic analysis code based on the small 
disturbance transonic equation and had only seven design 
variables. A more realistic problem would utilize an 
aerodynamic analysis code based on the full potential 
equation and might have fifteeen or more design variables. 
In this case the c.p.u. requirements of numerical 
optimization would generally be considered too large for the 
method to be used. 

It is found that a significant fraction of the time 
spent in designing airfoils using numerical optimization is 

spent on calculating the finite difference approximations to 
the gradients of the objective and constraint functions. 
Recall that a gradient calculation is made by separately 


perturbing each design 

variable 

and 

then 

performing 

an 

aerodynamic analysis 

to determine 

the 

change in 

the 

objective function and 

the active 

constraint 

functions. 

The 


flow about the unperturbed airfoil provides a good initial 
guess for the solution of the perturbed airfoil, but 
computer times are still very large. 

The objective of this research is to find a method for 
computing the gradients of the objective ana constraint 
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functions that require significantly less computer time than 
the method outlined above. The finite difference 
approximation to the gradients used only one sided 
differencing to approximate the gradients. One sided 
differencing is only first order accurate, suggesting that 
numerical optimization does not require extremely accurate 
gradient information. This research will take advantage of 
this fact and will try to make better use of the fact that 
the solution for the perturbed body is only slightly 
different from that of the unperturbed body. It is also an 
objective of this research to find a method to calculate 
gradients that is not specific to a particular set of 
governing equations. That is, the method should be 
applicable to the small disturbance transonic potential 
equation, the full potential equation, and Euler's 
equations . 
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CHAPTER 2. METHODS FOR FAST GRADIENT CALCULATIONS 

Five different methods were examined tha\ could quickly 
solve for the flow about an airfoil whose geometry has been 
slightly perturbed. One of these methods shows promise in 
meeting all the objectives stated above. The five methods 
are; 1 )Perturbation Equations, 2) Local Methods , 3) Local 
Linearization M) Method of Integral Relations, 5) Method of 
Reduced Domains. 

The first four methods listed were found to have various 
shortcomings while the final method listed shows some 
promise. This section will begin with a brief review of the 
first four methods listed above, after which a more 
extensive review of the Method of Reduced Domains will be 
given . 

All preliminary investigations were performed using the 
small disturbance transonic potential equation. A review of 
the important details involved in solving this equation is 
presented in the Appendix. 

1 )Perturbation Equations: 

An obvious first step in solving a flow problem which is 
a perturbation from a known base flow is to rewrite the 
governing equations with the parameters split into two 
parts. The first part would satisfy the base solution and 
the second part would be a perturbation from the base 
solution. Writing the potential as 

Cp-- 

-1 6 - 
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ana plugging this into the small disturbance potential 
equation yields 

Discarding higher order terms in the perturbation potential 
and noting that the base potential satisfies the equation 
yields the desired perturbation equation 

The perturbation form of the potential must also be 
substituted into the shock jump relations for a complete 
formulation of the problem. The method was not pursued, 
however, because equation 4 is essentially no simpler to 
solve than the original small disturbance equation. 


2) Local Methods 

Local methods try to relate the pressure at the airfoil 
surface with the local geometry of the airfoil. This is a 
method proposed by Davis (ref. 11) in which he uses the work 
of Spreiter and Alksne and their method of local 
linearization (ref. 12) as a basis to derive the following 
equations : 

Cp = ■ h nl M h / n, >1 (Sc) 


where 

nj : ('/(/- 
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and 


Cc 




I * ti* 



To obtain these very simple equations, a considerable 
number of questionable assumptions had to be made in 
addition to the original assumptions made by Spreiter and 
Alksne in their development of local linearization. Becuase 
of ail the assumptions made, the validity of equation 5 was 
put in doubt so the method was not pursued further. 


3) Local Linearization: 

An attempt was made to use the work done by Spreiter and 
Alksne that did not require the plethora of assumptions made 
in local methods. Local linearization is valid for purely 
subsonic flow, purely supersonic flow, and flow with free 
stream Mach number very close to 1. The method will be 
outlined for purely subsonic flow. The ideas used in this 
case are also usea in the supersonic case and the case where 
freestreani Mach number is near 1. 

The analysis begins with the small aisturbance transonic 
potential equation 

(i-n^-rC.Lfn)/p,)0,. > (Pyy ■ o U) 

Let 

'k = f?) 

and initially treat as a constant. This yielas the simple 
equation 
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the solution of which is 




This equation can now be differentiated with respect to x 


yielding 


J- Jt^L 


Oo) 


The expression for ^ is now substituted back into equation 
10 and this ordinary differential equation is solved 


yielding 


. 2 / 3 ^ (J^^C 


(-H) 


where 


k -■ ('>*1) 


and C is a constant of integration. The above step is an 
attempt to compensate for the approximation made initially 
in the analysis in which ^ is treated as a constant. The 
above steps may seem somewhat arbitrary, but Spreiter and 
Alksne have shown that this sequence of steps leads to the 
most reasonable approximation. 

The only thing that remains to be done is the evaluation of 
the constant of integration. The method given by Spreiter 
and Alksne for the evaluation of C is not used. Equation 11 
can be rewritten as 




The flow about a base airfoil is known so substituting in 
Ui and (/l in the above equation yielas an expression for 
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C. Note that C will now be a function of x. With C(x) 
known, u(x,0) can be found from equation 11 yielding 

Equation 13 was used to solve for Cp for the perturbed 
airfoil. The base airfoil used was a parabolic arc with S - 
.08. The freestreara Mach number was .75. This yielded a flow 
that was subsonic everywhere. Four perturbations in the 
form of shape functions were separately added to the base 
airfoil. The four shape functions are presented in figure 4. 
Each shape function was multiplied by a factor of .004. 

The flow about the perturbed airfoil was solved using 
equation 13 and also using the finite difference solver 
outlined in Appendix A. The changes in pressure coefficient 
from the base to perturbed airfoil for each perturbation is 
shown in figure 5. The results look very promising. The 
curves obtained using equation 13 or the finite difference 
method lie nearly on top of each other. The computer time 


required to 

solve the 

problem with equation 

13 

is 3 

orders 

of magnitude 
difference me 

less than 
thod . 

the time required 

by 

the 

finite 


Attempts were made to extend this method to 
supercritical flow, but these failed. The fundamental 
problem is that the ideas used in local linearization cannot 
be applied to mixed flow. The ability to extend this methoa 
to the full potential equation or Euler's equations also 
seems doubtful. 
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For these reasons, the method was abandonned. If, 
however, there is a need to design subcritical airfoils 
where the small aisturbance potential equation is a 
reasonable approximation, then this method should work very 
well . 


4) Method of Integral Relations 


The purpose of this method is to reduce the dimension of 

an equation by 1. How this is done is best explained by 

illustrating a simple example. Consider the following two 

dimensional equation 

^ f i o 

0 'y 

where the boundary conditions 


are given. This equation can be 


yielding 


Vi, 


*y 


^ 9* '?• = 

If F can be written as a(x)*b(x 
function then the integral can 


integrated with respect to y 

: o (/r) 

,y) where b(x,y) is a known 
be performed yielding 


(Bex) t ^ ^ 

0% 


where 



This one dimensional equation is much simpler to solve 
than the original two dimensional equation. The success of 
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this method depends on the accuracy with which b(x,y) is 
known. It was hoped that since a solution of a flow that is 
Just slightly perturbed from a base flow is desired, the 
base flow should provide a reasonable approximation to the 
function b(x,y). 

To use this method on the small oisturbance transonic 
potential equation, the equation must be written in 
divergence form as follows 

0^0.- i') ' ■■ ) 

thus in this case 

F- /<■ 4’A (/^‘■) 

G = <Py Cni.) 

0y is given at y = 0 and goes to 0 as 7 goes to infinity so 
the limits of integration are 0 and infinity. Writing F as 

gives for the function b(x,Y) 

^ } /fti ('<j o) 

The problem with this is that F(,(x,0) may be equal to zero 
at some point in the flow. In fact F^(x,7) way be zero at 
some point in the flow so that with b(x,T) written as a 
ratio, the possibility always exists that the denominator 
will be zero. 

Other forms for b(x,7) were investigated such as 
assuming that b(x,y) was an exponential function or an 
algebraic function. Again, the flow about the base airfoil 
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was used as a suide to determine the exact form that b(x,y) 
should take. No reasonable way of using the base airfoil to 
determine b(x,y) was found. For this reason the method coula 
not be used. 

5) Method of Reduced Domains 

The basic idea behind this method is very simple. When 
calculating the flow about the perturbed airfoil the same 
solver is used as the one used to calculate the flow about 
the base airfoil, but the size of the domain on which the 
calculation is made is greatly r ■educed. The motivation 
behind this is that a perturbation in airfoil geometry will 
primarily affect the flow very close to the airfoil so that 
the base flow will provide reasonable boundary conditions 
for the reduced domain. Figure 6 is a typical comparison of 
the domain size used to calculate the flow about the base 
airfoil r'ld the perturbed airfoil. 

The reduction in c r.puter time that can be expected 
using the method of reduced domains comes about not only 
t-.case the number of mesh points is fewer so the number of 
calculations per iteration is fewer, but also because fewer 
mesh points generally lead to a higher convi?rgence rate. 

Suppose that a problem is solved using an iterative, 
finite difference scheme and that we want to reduce the size 

* IFO 

cf the error by 10 . The number of iterations to ao this, 

p, is given by 

P- ^/Voq„0) 


(^ 1 ) 


. *✓ 
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where D Is the spectral radius. 
For a Laplace equation 

<P.. ' (Pyy ■- 




solved with a Jacobi iterative scheme on a square domain 
with Ax and Ay equal, we have 


9 = /- 




where N is the number of mesh points. 

The rate of convergence is thus dependent on the number 
of mesh points in this case. While the rate of convergence 
for the small disturbance potential equation or the full 
potential equation cannot be found analytically, it is 
expected that the rate of convergence will be a function of 
the number of mesh points. 

Unlike the first four methods tried this method makes 
its approximation in the boundary condtions used when 
solving for the flow about the perturbed airfoil, not in the 
actual method of solution. Some advantages to this are 
immediately apparent:!) Major changes to existing cooes are 
not required since the same solver is used, 2) Application 
of this method to codes based on different governing 
equations is possible, 3) A tradeoff exists between accuracy 
ana speed, the larger the oomain the more accurate the 
solution, but also the greater the c.p.u. requirement. This 
last advantage is important because it gives f le:<i'oility in 
using the method of reduced domains. In some applications 
the size of the ocmain coulc be reouced substatiaily thereby 
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greatly reducing c.p.u. requirements. Hov/ever, the size of 
the aomain could always be increased if accuracy became a 
problem. 

Some trends can be found in the accuracy of using the 
base flow for boundary conditions by examining the inner and 
outer expansions used at different Mach numbers to derive 
the small disturbance potential equation from the full 
potential equation. 

To derive the small disturbance potential equation valid 
for subsonic and supersonic flow 

(h * 0yy - O (l^) 

from the full potential equation 

C0‘- -f;) ' (o'- 4>y') 4>yy ' i iy <p.y - O U 


requi res 


the use of inner and outer expansions of the form 

I <k‘c^,y ) ' ^,7*. ?) * f,’ 4>.' (\ 9)*-1 


where 

y ■■ yA, 


and 6, is a small parameter related to the airfoil 
thickness . 

To derive the small disturbance potential equation valid 
for transonic flow 

[ / - ■ /?: (n<) (I'.] * Pyy -- O a?) 
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from the full potential equation requires the use of the 
following inner and outer expansions 

«« [ 4>> (>■ ')) ' fa <P:(T<.i) f (1 <Px <><■))' ■ ] 

where 

and f’i is a small parameter related to the airfoil 
thickness . 

The important point is that for supersonic and subsonic 
flow the inner solution is valid only for y of order €, , 
while for transonic flow the solution is valid for y of 
order 1. This means that a perturbation in airfoil geometry 
will affect the flow for a greater distance from the airfoil 
in transonic flow than in subsonic or supersonic flow. The 
accuracy of using the base flow for boundary conditions at 
the edge of a given reduced domain should therefore be less 
in transonic flow than in subsonic or supersonic flow. 

A complete outline of the inner and outer expansion 
procedure can be found in reference 13. 

Wind Tunnel Analogy 

There Is a physical atialog to computing the flow about 
the perturbed airfoil on a reduced domain that offers 
insight into the problem. Determining the aerodynamic 
characteristics of an airfoil by performing wine tunnel 
tests is similar to the methou of reduced domains. In eacli 
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case the flow is affected by the outer boundary. 

The hypothetical wind tunnel in which the perturbec 
airfoil is tested is, however, special. The solution fcr the 
base airfoil is obtained on the full domain so that if the 
domain is reduced, the exact boundary conditions are 
availiable at the edge of the reduced domain. If the flow in 
the interior of the reduced domaii. is perturbed and then the 
problem is solved with the same airfoil as before, exactly 
the same answer as that found on the full domain will be 
computed. This means that the hypothetical windtunnel has 
been constructed in such a way as to give aerodynamic 
characteristics for the base airfoil which are the same as 
if the base airfoil had been tested in free air. Figure 7 
shows the base airfoil in the hypothetical windtunnel. 
Solving for the flow about a slightly perturbed airfoil is 
thus like testing the slightly perturbed airfoil in the 
hypothetical windtunnel. 
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CHAPTER 3. RESULTS OF STUDIES USING THE METHOD OF REDUCED 

DOMAINS 

Exploratory Studies 

Results presented in this section compare the change in 
Cp caused by perturbations in airfoil geometry as computed 
on the full domain and the reduced domain. The base airfoil 
in each case was a parabolic arc at zero angle of attack, 
and the four perturbations of figure M multiplied by a small 
factor were separately added to the base airfoil. 

It is difficult to quantify the performance of the 

method of reduced domains without doing an optimization run. 
A subjective assessment of the results will be used 

initially in this chapter to develop confidence in the 

method. Later in this chapter an actual optimization test 
will be performed. 

Initially all results presented will be from subcritical 

tests. Greater difficulty was anticipated for supercritical 

\ 

tests and these are discussed after the subcritical results 
are properly understood. 

In solving the flow on the reduced domain, using as the 
basis for the solution the small disturbance potential 

equation, two different types of boundary conditions can be 
specified. Either the potential can be specified at the edge 
of the domain (Dirichlet boundary condition) or the normal 
derivative of the potential, the transverse velocity, can be 
specified (Neumann bounaary condition). 


» r 
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Solving the flow about the perturbeo airfoil on the 
reduced domain and specifying the potential to be that of 
the base flow on the outer boundary is analogous to testing 
the perturbed airfoil in the hypothetical windtunnel with a 
freejet boundary condition at the tunnel walls. That is, the 
specification of the potential to be that of the base flow 
on the outer boundary forces the pressure to be that of the 
base flow on the upper boundary of the reduced domain. In a 
windtunnel with free jet boundaries, one expects the peak 
perturbation velocities to be underestimated thus 
underestimating the peak perturbation Cp. 

Figure 8 shows the results where the base potential was 

specified at the edge of the reduced domain. The flow 

conditions in this case were: 

freestream Mach number =0.7 
thickness/chord = 0.1 

perturbation multiplication factor = 0.001 

The size of the reduced domain is shown in figure 9. The 

number of mesh points in this reduced domain was 116 while 

it was 3965 for the full domain. The underprediction of the 

change in Cp is consistent with the above argueraents. The 

results are, however , reasonably accurate and the average 

c.p.u. requirement in computing the flow on the reduced 

domain was approximentaly 60 times smaller than the c.p.u. 

requirement for computing the flow on the full comain. 

The next bounoary condition tested was a Neumann 
boundary condition on the upper boundary of the reduced 
domain. The finite difference algorithm used in these 
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computations required the specification of the potential at 
the upstream boundary, so the potential at the upstream 
boundary was specified to be that of the base flow. For 
convenience, the potential was specified to be that of the 
base flow at the downstream boundary also. In all 
computations where a Neumann boundary condition was used at 
the upper boundary, the above outlined boundary conditions 
were used at the upstream and downstream bounoaries. 

Specifying the transverse velocity to be that of the 
base flow on the edge of the reduced domain is analogous to 
solid wall boundaries on the hypothetical windtunnel. Solid 
wail boundaries on a windtunnel constrain the flow to be 
tangential to the solid wall. In this case one would expect 
a test in the hypothetical wind tunnel on the perturbed 
airfoil to overestimate peak velocities thus overestimating 
the perturbation in Cp. 

The results presented in figure 10 are consistent with 

this prediction. The flow conditions in this case were: 

freestream Mach number = 0.7 
thickness/chord = 0.1 

perturbation multiplication factor = 0.001 

The size of the domain used in this case is presented in 

figure 9. The average c.p.u. requirement for computing the 

flow on the reduced domain was reduced by a factor of 45 as 

compared to the c.p.u. requirement of the full domain. The 

smaller reduction in c.p.u. time in this case as compared to 

the previous case was anticipated since Neumann bounaary 

conaitions generally lead to slower convergence. 
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Windtunnels usually have ventilated walls to try to 
lessen the effect of the walls on the flow. The next section 
describes attempts to correct the boundary conditions used 
on the edge of the reduced domain in the hope of yieldig a 
more accurate solution. 

Boundary Condition Modifications 

Modifications were made to each of the boundary 
conditions used above. The Dirichlet boundary condition was 
modified by treating the perturbation in geometry of the 
airfoil as a perturbation doublet and then adding the 
perturbation in potential caused by this doublet to the base 
potential on the outer boundary. The Neumann boundary 
condition was modified by treating the perturbation in 
geometry of the airfoil as a wavy wall and using a 
simplified analysis to determine the effect of this 
perturbation on the transverse velocity at the boundary. 

In Appendix A it is shown that in the far field the 
airfoil is treated as though it were a doublet. This 
treatment of the airfoil allows the potential to be 
determined in the c'ar field by use of equation A5. This 
equation is accurate only if it is used at points a 
considerable distance from the airfoil ( at least 1 chord 
length ).The accuracy of equation A5 diminishes when it is 
applied closer and closer to the airfoil, but how rapidly 
the accuracy decays is not known. As an approximation, this 
equation was used to determine the change in potential at 
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the edge of the reduced domain caused by the perturbation in 
airfoil geometry. 

The doublet was positioned on the airfoil at the chord 
station where the amplitude of the perturbation was at a 
maximum. Equation A6 is the equation for the doublet 
strength. The equation has one part due to the airfoil 
volume and a second nonlinear part. The contribution from 
the nonlinear part is generally small so in computing the 
perturbation doublet strength it was negleted. The equation 


for the doublet strength was therefore 




( 19 ) 


where f (fj is the perturbation in airfoil geometery. 

Figure 11 presents the results obtained with the above 
outlined boundary conditions. The flow conditions for this 
test were: 

freestream Mach number =0.7 
thickness/chord = 0.1 

perturbation multiplication factor = 0.001 

The size of the reduced domain is given in figure 9. The two 

curves lie reasonably close to each other. The average 

computing time required on the reduced domain was 

approximately 40 times less than that required on the full 

domain . 


Modifications of Neumann boundary conditions is similar 
to changing the hypothetical windtunnel shape. The problem 
is to determine the extent of the modifications required. 
This problem can be restated as the problem of determining 
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how much of the perturbation in transverse velocity at the 
airfoil surface is translated to the edge of the reduced 
domain. As a guideline to determine the decay in the 
transverse velocity a simple wavy wall model for the 
perturbation is used. 

Consider the small disturbance transonic potential 
equation 

f - 0 C3o) 


subject to the boundary conditions 

(P. ; <Py are af oO 

and 


Let us solve this equation for subsonic flow where 
(/<’-( > ) is always greater than zero. To make the 

mathematics tractable let us approximate the coefficient of 
the 4>yy term as a constant, . This equation, with the 
boundary conditions given, can be solved in closed form 
using separation of variables. The solution is 

S> ( Cos C3/] 


Note that the local value of ^ is used at every point. 

The important information from equation 31 is the rate 
at which velocities at Y = 0 decay. Using the decay rate of 
this equation gives for the velocity at any point Y 




3i) 


- 35 - 



I 

I 



ORIGINAL {s 
OF POOR Quality 

Equation 32 can be used to approximate the perturbation 
velocity at the edge of the reduced domain. To use equation 
32 t (9 , which determines the wavelength of the wavy wall 
must be known. Since the airfoil is similar to a half sine 
wave, 6 was made to be TT . 

Figure 12 presents the results of using the 

approximations just oulined. The flow conditions in this 
test were: 

freestream Mach number = 0.7 

thickness/chord s 0.1 

perturbation multiplication factor = 0.001 

The size of the reduced domain is given by figure 9. The 
two curves lie very nearly on top of each other, indicating 
that the above approximation is a very good one. The average 
computing time required on the reduced domain was 50 times 
less than that required on the full domain. 

These subcritical results are encouraging. The results 
obtained with the four different boundary conditions used 
above appear to be at least reasonably accurate and 
sometimes very accurate. The saving in c.p.u. time in each 
case was substantial. 

The same four boundary conditions outlined above were 
tested in supercritical flow. The one difference is that the 
wavy wall formula given above is valid only in subsonic 
flow. A supersonic wavy wall formula is required for cases 
wnere the flow at the eage of the outer bounoary is 
supersonic. The derivation of this formula is completely 
analogous to the aerivation presented for the subsonic flow. 
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The result is 

As can be seen, the transverse velocity at the airfoil 
surface does not decay but propogates along characteristics. 
In locally supersonic flow the boundary condition was 
therefore applied as follows: The slope of the 

characteristic was determined from the local Mach number at 
the boundary. The characteristic was then approximated as a 
straight line. This line was followed to the airfoil surface 
where the perturbation transverse velocity was known. This 
transverse velocity was then propogated back up along the 
approximation of the characteristic to the boundary of the 
reduced domain, where it wa" added to the base transverse 
velocity. This procedure is illustrated in figure 13. 

The runs made with the boundary conditions being: 1) The 

base potential, 2) The base transverse velocity, and 3) The 

base potential with the potential due to a perturbation 

doublet added to it, were all rui' under the same conditions: 

freeesteam Mach number s 0.82 
thickness/chord = 0.1 

perturbation multiplication factor = 0.001 

The size of the domain used is shown in figure 9. 

Figure shows the results where the base potential was 
specified on the edge of the reduced domain. The effect of 
the perturbation doublet was negligible , the results being 
essentially those presented in figure 14. The results in 
these cases are not very gooa. The largest changes in Cp 
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I occur near the shock and these are not captured well. When 

the base transverse velocity was specified at the upper 
boundary of the reduced domain, the computations would not 

I 

converge to a solution. | 

"a 

The tests in subcritical flow using the wavy wall | 

approximation yielded the best results so it was hoped that 
they would yield good results in the supercritical case. The 
initial test was a conservative one in which the flow was 
only slightly supercritical. The test conditions in this 
case were: 

freestream Mach number = 0.78 
thickness/chord = 0.10 

perturbation multiplication factor = 0.0014 

The size of the reduced domain used is shown in figure 9. 

Results are presented in figure 15. As can be seen the 

results are not very accurate. It seems that supercritical 

flow is much more sensitive to boundary conditions than 

subcritical flow. 

The flexibility of the reduced domain method was put to 
use to try to overcome these difficulties. The outer edge of 
the domains used generally crossed through a region with 
supersonic flow. This was the probable cause of 

dissappointing results just presented. To get around this, 
the shape of the domain was altered so that the outer 
boundary was always in subsonic flow. A typical domain for 
supercritical flow calculations is presented in figure 16. 

Figure 16 also shows the boundary conditions tested. The 
base potential was specified everywhere except at the very 
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top of the "notch” part of the domain. There, either the 
base potential was specified or the base transverse velocity 
was specified. 


For each 

boundary 

condi tion 

a series of 

tests 

were 

conducted in 

which 

the domain 

size was increased 

until 

consistently 

acceptable results 

were obtained. 

For 

these 


tests the flow conditions were: 

freestream Mach number = 0.62 

thickness/chord s 0.10 

perturbation multipl icatin factor = 0.001 

Figure 17 presents what are considered acceptable 
results for the Dirichlet boundary condition and figure 18 
shows results for the Neumann boundary condition. The 
corresponding domain size required to achieve these results 
is shown in figure 19. The number of mesh points in tne 
reduced domain with Neumann boundary conditions was 216 
while the number of mesh points in the reduced domain with 

Dirichlet boundary conditions was ll^tO. Recall that the the 

number of mesh points in the full domain was 3965. As can 
be seen, the size of the domain required to obtain 
reasonable results using the Dirichlet boundary condition is 
considerably larger than the domain required to achieve 
reasonable results using the Neumann boundary condition. An 
explanation for this phenomena has not yet been found. 

The average c.p.u. time required to get the results on 
the reduced domain is about 30 times less chan on the full 
domain if a Neumann boundary condition is used and about 6 
times less if a Dirichlet boundary condition is used. While 
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the reduction in c.p.u. time is not as large as the 
reduction realized in subcritical flow, it is still 
substantial . 

Summary of Results of Exploratory Studies 

The results can be summarized as follows: 

Sucbritical Flow 

1) Reasonably accurate results can be obtained 
on a substantially reduced domain. 

2) There is flexibility in applying boundary 
conditions. Different boundary conditions 
yield different results but all of these 
results seem acceptable. 

3) The c.p.u. requirement of the reduced 
domain calculations is from 40 to 60 times 
lower than the c.p.u. requirement for full 
domain calculations, depending on boundary 
conditions applied. 

Supercritical Flow 

1) Reasonably accurate results can be obtained 
on a substantially reduced domain but the 
outer t undary of the reduced domain must be 
in subsonic flow. 

2) There is no flexibility in applying 
boundary conaitions. Fast and reliable results 
can only be obtained with the use of Neumann 
boundary conditions. 



3) The c. p. u. requirement of the reduced aoraain 
calculations is 30 times smaller than the 
c.p.u requirement of full oomain calculations. 

Note that the savings in computer time using the method 
of reduced domains came about primarily because of the 
reduction in the number of calculations required per 
iteration. The reduction in the actual number of iterations 
on the reduced domain was net large, always being less than 
a factor of 2. 

Full Potential Equation Calculations 

The knowledge gained from applying the method of reduced 
domains to an aerodynamic analysis code based on the small 
disturbance potential equation is now applied to the full 
potential equation. The full potential equation is 
considerably more accurate than the small disturbance 
potential equation and so is used much more frequently in 
the design of airfoils. 

A computer code that solves the full potential equation 
is much more complex than one that solves the small 
disturbance potential equation. FL06 , a nonconservative full 
potential equation solver written by Antony Jameson, was 
used in these tests. FL06 was chosen because it is the 
aerodynamic analysis code in use at NASA Ames for numerical 
optimization of airfoils. A complete outline of the 
fundamentals involved in solving the full potential equation 
can be found in reference Id. 




Testa conduct ed usinii FLOb tn which the size of the 
computational domain was greatly reduced. FLOb maps tlie 
airfoil onto a unit circle and then performs an Inversion 
such that Infinity Is mapped Into the origin of the unit 
circle. Thus the computational domain used In FL06 Is the 
Inside of a unit circle . The outer boundary of the reduced 
domain must be In subsonic flow, so for a typical case of a 
lifting airfoil wltij a supersonic zone on the upper surface 
of ti>e airfoil, the computational domain woula be like that 
siiown In figure JO. Because of the coordinate system used In 
FLOb, a reduction In domain size like tiie reduction shown In 
figure JO, reduced the number of mesh points by a factor of 
only 5 or b. 

FLOb also uses two steps 1 ti solving the equation: a 
relaxation step, and a fast solver step. The relaxation step 
gives slow convergence, but can be used on irregular domain 
shapes and Is valid for . upersonlc flow. The fast solver 
step gives very rapid convergence in subsonic flow, but 
becomes less and less stable as the local Mach number 
increases and goes unstable in supersonic, flow. Also, the 
fast solver can only be used on regular domain shapes. Tlie 
two steps combined converge reasonably fast In transonic 
f 1 ow . 

Because liio I'educed comput a 1 1 ona 1 doi.ialn is irr'eguiar in 
shape, the I’asl solver catinot be used over the wiiole domain. 
Therefore, the fast solver is used only in the region sliown 
In figiu'e JO. Tl'.e relaxtlon step is used over the whole 
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the fast solver 


on the 


domain. Not being able to uw.e 
complete reduced domain does not slow converge 
part of the domain where it is not used 
supersonic flow. 

Tests were run on the airfoil shown in fig 
angle of attack, with a freestream Mach n 
resulting in a supersonic zone on the u 
surface of the airfoil. A single perturbation 

Perl ' O.Ool * Sih(7J v) 


nee since that 
is a region of 

ure 21 at zero 
umber of 0.8 
pper and lower 
given by 

[if) 


was added to the upper and lower surface of the airfoil. The 
boundary conditions used at the edge of the reduced domain 
were: 1) Potential specified to be that of the base flow 
everywhere, 2) Potential specified to be that of the base 
flow everywhere except on the "notch” part of the grid where 
the normal velocity was specified to be that of the base 
flow. An unexpected result was found from these tests. While 
in case 1 the solution did not capture the shock movement as 
expected from tests performed on the small disturbance 
potential solver, case 2 converged so slowly that the 
required c.p.u. time was as large as that needed to solve 
for the perturbed flow on the full domain. 

This fundamental difference in results from the small 


disturbance 

potential solver and FL06 

did 

not 

seem 

reasonabl e . 

One important difference 

in the 

two 

codes 

was 

that FL06 

is a nonconservative 

code while 

the 

small 

disturbance 

solver is a conservative 

code. A 

compa 

rison 

was 
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made on the rate of convergence of the conservative ana 
nonconservative forms of the small disturbance code. These 
tests showed that while the two forms of the code had nearly 
identical convergence rates on the full domain, the 
conservative form of the code converged in fewer than half 
the number of iterations required by the nonconservative 
form of the code on the reduced domain. While not a 
conclusive result, this strongly suggests that the method of 
reduced domains can be made to work for supercritical flow 
if the aerodynamic analysis code is conservative. 

Optimization Calculation 

Since the method of reduced domains could not be made to 
work in supercritical flow using FL06, the method was tried 
on a subcritical case. The objective of this test was not to 
design a better airfoil but was to demonstrate that the 
method of reduced domains will produce reasonable results in 
an actual optimization test. 

The optimization code used for this test is called 
QNMDIFF. It employs a quasi -Newton method for unconstrained 
optimization. The base airfoil used for this test was a 17 
percent thick airfoil designed for general aviation 
applications (GA(W)-I). The airfoil geometry can be found in 
reference 15. Computations were made at a Mach number of 
0.22 and at zero angle of attack. The objective function 
used was the sum of the differences between a target 
pressure distribution along the upper surface of the airfoil 
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and the actual pressure distribut:* on, that is 

OBJ-- t fsr; 

The target pressure dlstrabutlon used was that of the actual 

(GA(W)-I) airfoil. The initial airfoil used was the 

(GA(W)-I) airfoil to which four shape functions had been 

added. Thus the optimized airfoil would be one in which the 

magnitude of the four shape functions was 0. The four shape 

functions used were the ones shown in figure 4. The initial 

magnitudes of the shape functions were: 

function 1 magnitude = 0.0100 
function 2 magnitude s 0.008? 
function 3 magnitude s 0.0093 
function 4 magnitude = 0.0021 

Table 1 compares the design vector and value of the 
objective function at the end of each optimization iteration 
as computed in the regular way and using the method of 
reduced domains. Figure 22 presents graphically the 

objective function history of the two methods. In each case 
the objective function had been reduced substantially after 
6 iterations and the design vector was very close to the 
exact answer of 0,0, 0,0. 

This optimization test required 40 aerodynamic 
evaluations in gradient calculations and 21 aerodynamic 
evaluations in line search calculations. QNMDIFF sometimes 
requires central differences in gradient calculations which 
is why the number of aerodynamic evaluations for gradient 
calculations is not 4»6=24. Each gradient calculation took 


time on the reduced 




approximately one tenth the c.p.u. 
domain as comparea to the full ooniain. The total amount of 
computer time required for this optimization calculation 
using the method of reduced domains was aproximately 45% of 
the time required in using the standard method. 
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The method of reduced domains has the potential to 
significantly reduce the computer time required in designing 
airfoils using numerical optimization. There are, however, 
unresolved complications in using the method if the flow is 
supercritical. The use of a nonconservative potential code 
leads to unacceptably slow convergence rates when 
computations are performed on the reduced domain. It is 
probable that this problem would be eliminated if a 
conservative potential code was used. 

Future work on the method of reduced domains should 
concentrate on understanding the effect that different 
boundary conditions have on computations using a reduced 
domain. This paper has presented results of using different 
boundary conditions. Explaining why one boundary condition 
works better than another has proved elusive. With the 
proper understanding of the effect of boundary conditions, 
the full potential of this method would be known. 

This research has concentrated on reducing the computer 
time required to determine the gradients needed by the 
optimization program. Research is also warranted on 
increasing the efficiency of the line search. Making full 
use of the fact that the line search is a one dimensional 
problem could yield a signifigant reauction in c.p.u. 
requirements. Another area of possible research is in 
determining the best way to modify the airfoils shape. 
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Different design variables will yield different design 
spaces and one design space may be much more conducive to 
numerical optimization than another. 
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Table 1. - Details of Optimization Test 











Figure 2. - Design Spec 
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Figure 18. - Perturbation in Cp as Computed on Full Domain and 
Reduced Domain. Base Transverse Velocity Specified 
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Figure 22. - Iteration History of Objective Function 
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APPENDIX. METHOD OF SOLUTION OF THE SMALL DISTURBANCE 
TRANSONIC POTENTIAL EQUATION 


Consider the similarity form of the small disturbance 
transonic potential equation 




y 


where 




ana 


y -■ 

The pressure coefficient is ther given by 

Cp - -i(p,(S Jr ) 




If the profile of the airfoil is given by 

y ' J M'i) 

then the Neumann boundary condition, transi'errec to the axis 

is . . ^ 

(P ^ ^ J 

In the far field the airfoii is treated as a doublet. This 
gives for the potential in the far field 

where () = doublet stength ^ 

: 3 jJ(J) >■ If dj (Ai) 


«o 
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The numerical method used to solve this equation 
proceeds as follows. Let p be a central difference 
approimation to the x derivatives 


where 


^ ‘i u [ ~ ^ j 

^ J 



and let be a central difference approximation to the y 
derivatives 

C/U) 

■ay' 

Define the following switching function 

z O if /l-'.j >0 (sUso^,^) (A 9a) 

Jti,^ - / >f (A9h) 

Then a fully conservative finite difference scheme to solve 
the small disturbance potential equa is 


Details of this method can be found in references 14,16, 


and 17. 
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